i = 0
nn.results <- list()
size_seq <- seq(4,12,2)
  i = i + 1
  # Fit GBM
  nn.results[[i]] <- list()
  N <- length(dta.train[,dv])
  cv_results <- foreach (cv=1:cv.runs, 
                         .combine='rbind',
                         .packages = c("nnet",
                                       "pROC")) %dopar% {   # Loop over CV runs
                 set.seed(52*cv) 
                 shuffle <- sample(1:N,
                                   N,
                                   replace=FALSE)
                 pred.prob <- matrix(NA,nrow=N,ncol=length(size_seq))
                 sampleSplit = c(0,round((N/5*c(1:5))))
                 for (f in 1:5) {
                   # Remove constant variables from predictors
                   train.RHS <- dta.train[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                 rhs]
                   test.RHS <- dta.train[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                              rhs]
                   not.constant <- apply(train.RHS,2,sd)!=0
                   train.RHS <- train.RHS[,not.constant]
                   test.RHS <- test.RHS[,not.constant]
                   # Find PC weights
                   pcs.fit <- prcomp(train.RHS,center = TRUE, scale = TRUE)
                   train.pcs <- predict(pcs.fit,
                                        train.RHS)[,1:min(nn_params$pcNum[[country]],
                                                          ncol(train.RHS))]
                   test.pcs <- predict(pcs.fit,
                                       test.RHS)[,1:min(nn_params$pcNum[[country]],
                                                        ncol(train.RHS))]
                   for (s in 1:length(size_seq)) {
                     # Fit Neural Network
                     cv.fit = nnet(x = train.pcs,
                                   y = as.matrix(dta.train[-shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                                                               dv]),
                                   size = size_seq[s], 
                                   decay = nn_params$decay,
                                   trace = nn_params$trace
                     )
                     pred.prob[shuffle[(sampleSplit[f]+1):(sampleSplit[f+1])],
                               s] <- predict(cv.fit,
                                            newdata = as.matrix(test.pcs))
                   }
                 }
                 cv_results <- c("CV Run"=0,
                                 "Size"=0,
                                 "AUC"=0)
                 # Get error rates for each lambda used
                 for (s in 1:length(size_seq)) {
                     auc <- roc(as.matrix(dta.train[,dv]),
                                pred.prob[,s])$auc
                     if (auc>cv_results["AUC"]) {
                       cv_results <- c("CV Run"=cv,
                                       "Size" = size_seq[s],
                                       "AUC" = auc)
                       best.s <- s
                     }
                     print(auc)
                 }
                 cv_results<- c(cv_results,
                                "DT" = as.numeric(mostaccDT(y=dta.train[,dv],
                                                            yhat=pred.prob[,s],
                                                            J=1000)))
                 return(cv_results)
  }
  hours <- floor((proc.time()[3]-start.time)/3600)
  mins <- floor((proc.time()[3]-start.time)/60) - hours*60
  secs <- floor((proc.time()[3]-start.time)) - hours*3600 - mins*60
  
  print(paste(hours,
              "h",
              mins,
              "m",
              secs,
              "s elapsed",
              sep = ""))
  
  nn.results[[i]]$size <- round(median(cv_results[,"Size"]))
  nn.results[[i]]$dt <- mean(cv_results[,"DT"])
  dta.train.rhs <- dta.train[,rhs]
  not.constant <- apply(dta.train.rhs,2,sd)!=0
  dta.train.rhs <- dta.train.rhs[,not.constant]
  pcs.fit <- prcomp(dta.train.rhs,
                    center = TRUE, 
                    scale = TRUE)
  insample.pcs <- predict(pcs.fit,dta.train.rhs)[,1:min(nn_params$pcNum[[country]],ncol(dta.train.rhs))]
  outsample.pcs <- predict(pcs.fit,dta.test)[,1:min(nn_params$pcNum[[country]],ncol(dta.train.rhs))]
  
  set.seed(i)
  nn.results[[i]]$mod <- nnet(x = insample.pcs,
                   y = as.matrix(dta.train[,
                                          dv]),
                   size = nn.results[[i]]$size,
                   decay = nn_params$decay,
                   trace = nn_params$trace
  )
  nn.results[[i]]$fit.oos <- predict(nn.results[[i]]$mod,
                                      newdata = outsample.pcs)
  nn.results[[i]]$actual.oos <- as.matrix(dta.test[,
                                                    dv])

  save(nn.results,file=paste(modeldir,"/",
                             country,
                             "_nn_",
                             v,
                             "_",
                             rhs.group,
                             "_cross",
                             fileext,
                             ".RData",
                             sep = ""))